C
      SUBROUTINE BTN4DF(INBTN,IOUT,ISUM,ISUM2,NCOL,NROW,NLAY,NPER,
     & NCOMP,MCOMP,TRNOP,TUNIT,LUNIT,MUNIT,NODES,MXCOMP)
C ****************************************************************
C THIS SUBROUTINE READS PROBLEM DIMENSIONS AND TRANSPORT OPTIONS.
C ****************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   INBTN,IOUT,ISUM,ISUM2,NCOL,NROW,NLAY,NPER,
     &          NCOMP,MCOMP,I,NODES,IERR,MXCOMP
      LOGICAL   TRNOP(10)
      CHARACTER HEADNG(2)*80,TUNIT*4,LUNIT*4,MUNIT*4
C
C--READ AND PRINT HEADING
      READ(INBTN,'(A80)') (HEADNG(I),I=1,2)
      WRITE(IOUT,10)
      WRITE(IOUT,15) (HEADNG(I),I=1,2)
      WRITE(IOUT,10)
   10 FORMAT(1X,' ----- ')
   15 FORMAT(1X,'| M T | ',A80/1X,'| 3 D | ',A80)
C
C--READ AND PRINT NO. OF LAYERS, ROWS, COLUMNS, AND STRESS PERIODS,
C--COMPONENTS
      READ(INBTN,'(6I10)',ERR=25,IOSTAT=IERR)
     &                    NLAY,NROW,NCOL,NPER,NCOMP,MCOMP
      IF(NCOMP.LT.1) NCOMP=1
      IF(MCOMP.LT.1) MCOMP=1
      IF(NCOMP.GT.MXCOMP.OR.MCOMP.GT.MXCOMP) THEN
        WRITE(*,1005)
		call stopfile  ! emrl
        STOP
      ENDIF
   25 IF(IERR.NE.0) THEN
        BACKSPACE(INBTN)
        READ(INBTN,'(4I10)') NLAY,NROW,NCOL,NPER
        NCOMP=1
        MCOMP=1
      ENDIF
      WRITE(IOUT,1000) NLAY,NROW,NCOL,NPER,NCOMP,MCOMP
 1000 FORMAT(1X,'THE TRANSPORT MODEL CONSISTS OF ',I5,' LAYER(S)',I5,
     & ' ROW(S)',I5,' COLUMN(S)',
     & /1X,'NUMBER OF STRESS PERIOD(S) FOR TRANSPORT SIMULATION =',I5,
     & /1X,'NUMBER OF ALL COMPONENTS INCLUDED IN SIMULATION =',I5,
     & /1X,'NUMBER OF MOBILE COMPONENTS INCLUDED IN SIMULATION =',I5)
 1005 FORMAT(/1X,'ERROR: MAXIMUM NUMBER OF COMPONENTS EXCEEDED!',
     &       /1X,'INCREASE DIMENSION OF [MXCOMP] IN THE MAIN PROGRAM.')
C
C--READ AND PRINT UNITS FOR TIME, LENGTH AND MASS TO BE USED
      READ(INBTN,'(3A4)') TUNIT,LUNIT,MUNIT
      WRITE(IOUT,1010) TUNIT,LUNIT,MUNIT
 1010 FORMAT(1X,'UNIT FOR TIME IS ',A4,';',2X,'UNIT FOR LENGTH IS ',
     & A4,';',2X,'UNIT FOR MASS IS ',A4)
C
C--READ AND PRINT MAJOR TRANSPORT OPTIONS USED IN THE SIMULATION
      DO I=1,10
        TRNOP(I)=.FALSE.
      ENDDO
      READ(INBTN,'(10L2)',ERR=100) (TRNOP(I),I=1,10)
  100 WRITE(IOUT,1020)
      WRITE(IOUT,'(1X,10I3)') (I,I=1,10)
      WRITE(IOUT,'(1X,10L3)') (TRNOP(I),I=1,10)
      WRITE(IOUT,'(1X)')
 1020 FORMAT(1X,'PACKAGES INCLUDED IN CURRENT SIMULATION:')
C
C--INITIALIZE ARRAY POINTERS FOR ALLOCATING MEMORY
      ISUM=1
      ISUM2=1
C
C--GET TOTAL NUMBER OF MODEL NODES
      NODES=NCOL*NROW*NLAY
C
C--NORMAL RETURN
      RETURN
      END
C
C
      SUBROUTINE BTN4AL(INBTN,IOUT,ISUM,ISUM2,NCOL,NROW,NLAY,NCOMP,
     & LCLAYC,LCDELR,LCDELC,LCHTOP,LCDZ,LCPR,LCXBC,LCYBC,LCZBC,
     & LCQX,LCQY,LCQZ,LCQSTO,LCDH,LCIB,LCCOLD,LCCNEW,LCCWGT,
     & LCCADV,LCRETA,LCSR,LCBUFF)
C **********************************************************************
C THIS SUBROUTINE ALLOCATES SPACE FOR ARRAYS NEEDED BY THE ENTIRE MODEL.
C **********************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   INBTN,IOUT,ISUM,ISUM2,NCOL,NROW,NLAY,NCOMP,LCSR,
     & LCLAYC,LCDELR,LCDELC,LCHTOP,LCDZ,LCPR,LCXBC,LCYBC,
     & LCZBC,LCQX,LCQY,LCQZ,LCDH,LCIB,LCCOLD,LCCNEW,LCCWGT,LCCADV,
     & LCRETA,LCBUFF,NODES,ISUMX,ISUMIX,ISOLD,ISOLD2,LCQSTO
C
C--PRINT PACKAGE NAME AND VERSION NUMBER
      WRITE(IOUT,1030) INBTN
 1030 FORMAT(1X,'BTN4 -- BASIC TRANSPORT PACKAGE,',
     & ' VERSION 4, AUGUST 2001, INPUT READ FROM UNIT',I3)
C
C--ALLOCATE SPACE FOR ARRAYS
      ISOLD=ISUM
      ISOLD2=ISUM2
      NODES=NCOL*NROW*NLAY
C
C--INTEGER ARRAYS
      LCLAYC=ISUM2
      ISUM2=ISUM2+NLAY
      LCIB=ISUM2
      ISUM2=ISUM2+NODES * NCOMP
C
C--REAL ARRAYS
      LCDELR=ISUM
      ISUM=ISUM+NCOL
      LCDELC=ISUM
      ISUM=ISUM+NROW
      LCHTOP=ISUM
      ISUM=ISUM+NCOL*NROW
      LCDZ=ISUM
      ISUM=ISUM+NODES
      LCPR=ISUM
      ISUM=ISUM+NODES
      LCXBC=ISUM
      ISUM=ISUM+NCOL
      LCYBC=ISUM
      ISUM=ISUM+NROW
      LCZBC=ISUM
      ISUM=ISUM+NODES
      LCQX=ISUM
      ISUM=ISUM+NODES
      LCQY=ISUM
      ISUM=ISUM+NODES
      LCQZ=ISUM
      ISUM=ISUM+NODES
      LCQSTO=ISUM
      ISUM=ISUM+NODES
      LCDH=ISUM
      ISUM=ISUM+NODES
      LCCOLD=ISUM
      ISUM=ISUM+NODES * NCOMP
      LCCNEW=ISUM
      ISUM=ISUM+NODES * NCOMP
      LCCWGT=ISUM
      ISUM=ISUM+NODES * NCOMP
      LCCADV=ISUM
      ISUM=ISUM+NODES * NCOMP
      LCRETA=ISUM
      ISUM=ISUM+NODES * NCOMP
      LCSR=ISUM
      ISUM=ISUM+NODES * NCOMP
      LCBUFF=ISUM
      ISUM=ISUM+NODES
C
C--CHECK HOW MANY ELEMENTS OF THE X AND IX ARRAYS ARE USED
      ISUMX=ISUM-ISOLD
      ISUMIX=ISUM2-ISOLD2
      WRITE(IOUT,1090) ISUMX,ISUMIX
 1090 FORMAT(1X,I10,' ELEMENTS OF THE  X ARRAY USED BY THE BTN PACKAGE'
     & /1X,I10,' ELEMENTS OF THE IX ARRAY USED BY THE BTN PACKAGE'/)
C
C--NORMAL RETURN
      RETURN
      END
C
C
C emrl - added ispc,species to subroutine below
      SUBROUTINE BTN4RP(IN,IOUT,IUCN,IOBS,IMAS,ICNF,ICBM,ispc,species,
     & NCOL,NROW,NLAY,NCOMP,LAYCON,DELR,DELC,HTOP,DZ,PRSITY,ICBUND,COLD,
     & CNEW,CADV,CINACT,THKMIN,XBC,YBC,ZBC,RETA,RFMIN,BUFF,MXPRS,NPRS,
     & TIMPRS,MXOBS,NOBS,NPROBS,LOCOBS,TUNIT,LUNIT,MUNIT)
C **********************************************************************
C THIS SUBROUTINE READS AND PREPARES INPUT DATA RELEVANT TO THE ENTIRE
C SIMULATION.
C***********************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   IN,IOUT,IUCN,IOBS,IMAS,ICNF,ICBM,NCOL,NROW,NLAY,LAYCON,
     &          NCOMP,ICBUND,MXPRS,NPRS,MXOBS,NOBS,LOCOBS,IFMTCN,IFMTNP,
     &          IFMTRF,IFMTDP,N,J,I,K,IP1,IERR,NPROBS,NPRMAS,INDEX
      REAL      DELR,DELC,HORIGN,HTOP,DZ,PRSITY,CNEW,COLD,CINACT,TIMPRS,
     &          XBC,YBC,ZBC,ZZ,XMAX,YMAX,ZMAX,RETA,RFMIN,BUFF,TEMP,CADV,
     &          THKMIN,CDRY
      LOGICAL   UNIFOR,UNIDX,UNIDY,UNIDZ,SAVUCN,SAVCBM,CHKMAS
      CHARACTER ANAME*24,FLNAME*40,FINDEX*30,TUNIT*4,LUNIT*4,MUNIT*4 ! emrl flname 50to40
      DIMENSION LAYCON(NLAY),ICBUND(NCOL,NROW,NLAY,NCOMP),DELR(NCOL),
     &          DELC(NROW),DZ(NCOL,NROW,NLAY),PRSITY(NCOL,NROW,NLAY),
     &          HTOP(NCOL,NROW),XBC(NCOL),YBC(NROW),ZBC(NCOL,NROW,NLAY),
     &          CNEW(NCOL,NROW,NLAY,NCOMP),COLD(NCOL,NROW,NLAY,NCOMP),
     &          RETA(NCOL,NROW,NLAY,NCOMP),BUFF(NCOL,NROW,NLAY),
     &          CADV(NCOL,NROW,NLAY,NCOMP),TIMPRS(MXPRS),LOCOBS(3,MXOBS)
      COMMON /PD/HORIGN,XMAX,YMAX,ZMAX,UNIDX,UNIDY,UNIDZ
      COMMON /OC/IFMTCN,IFMTNP,IFMTRF,IFMTDP,SAVUCN,SAVCBM,CHKMAS,NPRMAS
C-------emrl start
	character*256 ext,species(20),flname2  !emrl 80t256
	integer      ISPC,IVERSION,IOBTY,N6,ISFLT,JSFLT,ISFLG,JSFLG,ISCL,
     &             INODE,IELEM,NNP,INAME
C-------emrl end
C
C--READ INPUT DATA
C  ===============
C
C--READ AND ECHO LAYER TYPE CODES
      READ(IN,'(40I2)') (LAYCON(K),K=1,NLAY)
      WRITE(IOUT,1000)
 1000 FORMAT(1X,'LAYER NUMBER  AQUIFER TYPE',
     &      /1X,'------------  ------------')
      DO K=1,NLAY
        WRITE(IOUT,1010) K,LAYCON(K)
      ENDDO
 1010 FORMAT(1X,4X,I3,10X,I3)
C
C--CALL RARRAY TO READ IN CELL WIDTH ALONG ROWS
      ANAME='WIDTH ALONG ROWS (DELR)'
      CALL RARRAY(DELR(1),ANAME,1,NCOL,0,IN,IOUT)
C
C--CHECK WHETHER ELEMENTS OF DELR ARE UNIFROM
      UNIDX=UNIFOR(DELR(1),NCOL,1,1)
C
C--CALL RARRAY TO READ IN CELL WIDTH ALONG COLUMNS
      ANAME='WIDTH ALONG COLS (DELC)'
      CALL RARRAY(DELC(1),ANAME,1,NROW,0,IN,IOUT)
C
C--CHECK WHETHER ELEMENTS OF DELC ARE UNIFROM
      UNIDY=UNIFOR(DELC(1),NROW,1,1)
C
C--CALL RARRAY TO READ IN TOP ELEVATION OF 1ST LAYER
      ANAME='TOP ELEV. OF 1ST LAYER'
      CALL RARRAY(HTOP,ANAME,NROW,NCOL,0,IN,IOUT)
C
C--CALL RARRAY TO READ IN THICKNESS ONE LAYER AT A TIME
      ANAME='CELL THICKNESS (DZ)'
      DO K=1,NLAY
        CALL RARRAY(DZ(1,1,K),ANAME,NROW,NCOL,K,IN,IOUT)
      ENDDO
C
C--CHECK WHETHER VERTICAL DISCRTIZATION IS HORIZONTAL
      UNIDZ=UNIFOR(HTOP,NCOL,NROW,1)
     & .AND.UNIFOR(DZ,NCOL,NROW,NLAY)
C
C--CALL RARRAY TO READ IN POROSITY ONE LAYER AT A TIME
      ANAME='POROSITY'
      DO K=1,NLAY
        CALL RARRAY(PRSITY(1,1,K),ANAME,NROW,NCOL,K,IN,IOUT)
      ENDDO
C
C--CALL IARRAY TO READ IN CONCENTRATION BOUNDARY ARRAY
      ANAME='CONCN. BOUNDARY ARRAY'
      DO K=1,NLAY
        CALL IARRAY(ICBUND(1,1,K,1),ANAME,NROW,NCOL,K,IN,IOUT)
      ENDDO
C
C--CALL RARRAY TO READ IN INITIAL CONCENTRATION
      ANAME='INITIAL CONC.: COMP. NO.'
      DO INDEX=1,NCOMP
        WRITE(ANAME(22:24),'(I3.2)') INDEX
        DO K=1,NLAY
          CALL RARRAY(COLD(1,1,K,INDEX),ANAME,NROW,NCOL,K,IN,IOUT)
        ENDDO
      ENDDO
C
C--READ AND ECHO CINACT,THKMIN
      READ(IN,'(2F10.0)',ERR=50,IOSTAT=IERR) CINACT,THKMIN
      IF(THKMIN.LT.0) THKMIN=0.
   50 IF(IERR.NE.0) THEN
        BACKSPACE (IN)
        READ(IN,'(F10.0)') CINACT
        THKMIN=0.
      ENDIF
      WRITE(IOUT,1020) CINACT,THKMIN
      IF(THKMIN.GT.0.05) THEN
        WRITE(IOUT,1022)
        THKMIN=0.01
      ENDIF
 1020 FORMAT(/1X,'VALUE INDICATING INACTIVE CONCENTRATION CELLS = ',
     & G15.7/1X,'MINIMUM SATURATED THICKNESS [THKMIN] ',
     & 'ALLOWED =',F8.4,' OF TOTAL CELL THICKNESS')
 1022 FORMAT(1X,'[THKMIN] MUST BE < OR = 0.05;',
     & ' RESET TO DEFAULT VALUE OF 0.01 [1% OF TOTAL CELL THICKNESS]')
C
C--READ AND ECHO OUTPUT CONTROL OPTIONS
      READ(IN,'(4I10,L10)') IFMTCN,IFMTNP,IFMTRF,IFMTDP,SAVUCN
      SAVCBM=.FALSE.
      WRITE(IOUT,1025)
C
      IF(IFMTCN.NE.0) WRITE(IOUT,1030) IFMTCN
      IF(IFMTCN.EQ.0) WRITE(IOUT,1032)
 1025 FORMAT(//1X,'OUTPUT CONTROL OPTIONS'/1X,22('-'))
 1030 FORMAT(/1X,'PRINT CELL CONCENTRATION USING FORMAT CODE:',I5)
 1032 FORMAT(/1X,'DO NOT PRINT CELL CONCENTRATION')
C
      IF(IFMTNP.NE.0) WRITE(IOUT,1034) IFMTNP
      IF(IFMTNP.EQ.0) WRITE(IOUT,1036)
 1034 FORMAT(1X,'PRINT PARTICLE NUMBER IN EACH CELL',
     &          ' USING FORMAT CODE:',I5)
 1036 FORMAT(1X,'DO NOT PRINT PARTICLE NUMBER IN EACH CELL')
C
      IF(IFMTNP.NE.0) WRITE(IOUT,1038) IFMTRF
      IF(IFMTNP.EQ.0) WRITE(IOUT,1040)
 1038 FORMAT(1X,'PRINT RETARDATION FACTOR USING FORMAT CODE:',I5)
 1040 FORMAT(1X,'DO NOT PRINT RETARDATION FACTOR')
C
      IF(IFMTDP.NE.0) WRITE(IOUT,1042) IFMTDP
      IF(IFMTDP.EQ.0) WRITE(IOUT,1044)
 1042 FORMAT(1X,'PRINT DISPERSION COEFFICIENT USING FORMAT CODE:',I5)
 1044 FORMAT(1X,'DO NOT PRINT DISPERSION COEFFICIENT')
C
c emrl      IF(SAVUCN) THEN
c emrl        WRITE(IOUT,1046) IUCN
c emrl        FLNAME='MT3Dnnn.UCN'
c emrl        DO INDEX=1,NCOMP
c emrl          WRITE(FLNAME(5:7),'(I3.3)') INDEX
c emrl          CALL OPENFL(-(IUCN+INDEX),0,FLNAME,1,FINDEX)
c emrl        ENDDO
c emrl      ELSE
c emrl        WRITE(IOUT,1047)
c emrl      ENDIF
c-------emrl start
      IF(SAVUCN) THEN
          IVERSION = 3000
          IOBTY = 100
          N6 = 7
          ISFLT = 110
          JSFLT = 4
          ISFLG = 120
          JSFLG = 4
          ISCL = 130
          INODE = 170
          IELEM = 180
          NNP = NROW*NCOL*NLAY
          INAME = 190
        DO INDEX=1,NCOMP
          WRITE(ISPC+INDEX) IVERSION,IOBTY,N6,ISFLT,JSFLT,ISFLG,JSFLG
          WRITE(ISPC+INDEX) ISCL
          WRITE(ISPC+INDEX) INODE,NNP
          WRITE(ISPC+INDEX) IELEM,NNP
C-----------emrl - always make sure that flname is 40 char
          FLNAME = SPECIES(INDEX)
          WRITE(ISPC+INDEX) INAME,FLNAME
        ENDDO
      ELSE
        WRITE(IOUT,1047)
      ENDIF
c-------emrl end
 1046 FORMAT(1X,'SAVE CONCENTRATION IN UNFORMATTED FILE ',
     & '[MT3Dnnn.UCN] ON UNIT ',I3)
 1047 FORMAT(1X,'DO NOT SAVE CONCENTRATION IN UNFORMATTED FILE')
C
C--READ NUMBER OF TIMES AT WHICH SIMULATION RESULTS SHOULD BE SAVED
C--IN STANDARD OUTPUT FILE OR RECORDED IN UNFORMATTED FILE
      READ(IN,'(I10)') NPRS
      IF(NPRS.LT.0) THEN
        WRITE(IOUT,1050) -NPRS
      ELSE
        WRITE(IOUT,1052) NPRS
      ENDIF
      IF(NPRS.GT.MXPRS) THEN
        WRITE(*,1054) MXPRS
		call stopfile  ! emrl
        STOP
      ENDIF
      IF(NPRS.GT.0) THEN
        READ(IN,'(8F10.0)') (TIMPRS(I),I=1,NPRS)
C
C--MAKE SURE ELEMENTS IN ARRAY [TIMPRS] ARE MONOTONICALLY INCREASING
        DO I=1,NPRS-1
          DO IP1=I+1,NPRS
            IF(TIMPRS(I).GT.TIMPRS(IP1)) THEN
              TEMP=TIMPRS(I)
              TIMPRS(I)=TIMPRS(IP1)
              TIMPRS(IP1)=TEMP
            ENDIF
          ENDDO
        ENDDO
        WRITE(IOUT,1055) (TIMPRS(I),I=1,NPRS)
      ENDIF
 1050 FORMAT(/1X,'SIMULATION RESULTS ARE SAVED EVERY ',I3,
     & ' TRANSPORT STEP(S)')
 1052 FORMAT(/1X,'NUMBER OF TIMES AT WHICH SIMULATION RESULTS',
     &' ARE SAVED =',I5)
 1054 FORMAT(/1X,'ERROR: MAXIMUM NUMBER OF TIMES AT WHICH',
     & ' SIMULATION RESULTS CAN BE SAVED IS',I5,
     & /1X,'INCREASE DIMENSION OF [MXPRS] IN THE MAIN PROGRAM')
 1055 FORMAT(1X,'TOTAL ELAPSED TIMES AT WHICH SIMULATION RESULTS ',
     & 'ARE SAVED: ',100(/1X,8G13.5))
C
C--READ NUMBER OF OBSERVATION POINTS
      READ(IN,'(2I10)',ERR=100,IOSTAT=IERR) NOBS,NPROBS
      IF(NPROBS.LT.1) NPROBS=1
  100 IF(IERR.NE.0) THEN
        BACKSPACE (IN)
        READ(IN,'(I10)') NOBS
        NPROBS=1
      ENDIF
      WRITE(IOUT,1056) NOBS
      IF(NOBS.GT.MXOBS) THEN
        WRITE(*,1058) MXOBS
		call stopfile  ! emrl
        STOP
      ENDIF
      IF(NOBS.GT.0) THEN
        WRITE(IOUT,1062) IOBS,NPROBS
        WRITE(IOUT,1060)
        DO N=1,NOBS
          READ(IN,'(3I10)') (LOCOBS(I,N),I=1,3)
          WRITE(IOUT,'(I4,4X,3(I5,2X))') N,(LOCOBS(I,N),I=1,3)
        ENDDO
        FLNAME='MT3Dnnn.OBS'
        DO INDEX=1,NCOMP
          WRITE(FLNAME(5:7),'(I3.3)') INDEX
          CALL OPENFL(IOBS+INDEX,0,FLNAME,1,FINDEX)
          WRITE(IOBS+INDEX,1063) ((LOCOBS(I,N),I=1,3),N=1,NOBS)
        ENDDO
      ENDIF
 1056 FORMAT(/1X,'NUMBER OF OBSERVATION POINTS =',I5)
 1058 FORMAT(/1X,'ERROR: MAXIMUM NUMBER OF OBSERVATION POINTS IS',I5,
     & /1X,'INCREASE DIMENSION OF [MXOBS] IN THE MAIN PROGRAM')
 1060 FORMAT(1X,'LOCATION OF OBSERVATION POINTS'/1X,30('.')
     &      /1X,'NUMBER  LAYER   ROW   COLUMN')
 1062 FORMAT(1X,'CONCENTRATION AT OBSERVATION POINTS SAVED IN FILE',
     & ' [MT3Dnnn.OBS] ON UNIT ',I3,' EVERY',I4,' TRANSPORT STEPS')
 1063 FORMAT(1X,' STEP   TOTAL TIME',
     & '             LOCATION OF OBSERVATION POINTS (K,I,J)'
     & /1X,17X,16(1X,3I4,1X)/(1X,17X,16(1X,3I4,1X)))
C
C--READ AND ECHO LOGICAL FLAG CHKMAS
      READ(IN,'(L10,I10)',ERR=105,IOSTAT=IERR) CHKMAS,NPRMAS
      IF(NPRMAS.LT.1) NPRMAS=1
  105 IF(IERR.NE.0) THEN
        BACKSPACE (IN)
        READ(IN,'(L10)') CHKMAS
        NPRMAS=1
      ENDIF
      IF(CHKMAS) THEN
        WRITE(IOUT,1064) IMAS,NPRMAS
c emrl        FLNAME='MT3Dnnn.MAS'
        ext = 'mas'  ! emrl
        DO INDEX=1,NCOMP
	    call setname(ext,flname2,index)  ! emrl
c emrl          WRITE(FLNAME(5:7),'(I3.3)') INDEX
c emrl          CALL OPENFL(IMAS+INDEX,0,FLNAME,1,FINDEX)
          CALL OPENFL(IMAS+INDEX,0,flname2,1,FINDEX)
          WRITE(IMAS+INDEX,1066)
          WRITE(IMAS+INDEX,1068) TUNIT,MUNIT,MUNIT,MUNIT,MUNIT
        ENDDO
      ELSE
        WRITE(IOUT,1065)
      ENDIF
 1064 FORMAT(/1X,'A ONE-LINE SUMMRY OF MASS BALANCE SAVED IN FILE',
     & ' [MT3Dnnn.MAS] ON UNIT ',I3,' EVERY',I4,' TRANSPORT STEPS')
 1065 FORMAT(/1X,'A ONE-LINE SUMMRY OF MASS BALANCE',
     & ' WILL NOT BE SAVED')
 1066 FORMAT(1X,'     TIME       TOTAL IN      TOTAL OUT      SOURCES',
     & '        SINKS      NET MASS FROM  TOTAL MASS',
     & '        DISCREPANCY(%)')
 1068 FORMAT(1X,5(4X,'(',A4,')',4X),
     & ' FLUID-STORAGE  IN AQUIFER  (TOTAL IN-OUT)  (ALTERNATIVE)')
C
C--SAVE MODEL GRID CONFIGURATION IN FILE [MT3D.CNF]
C--FOR USE WITH UNFORMATTED CONCENTRATION FILE BY POST-PROCESSOR
c emrl      IF(SAVUCN) THEN
c emrl        CDRY=CINACT
c emrl        FLNAME='MT3D.CNF'
c emrl        CALL OPENFL(ICNF,0,FLNAME,1,FINDEX)
c emrl        WRITE(ICNF,*) NLAY,NROW,NCOL
c emrl        WRITE(ICNF,*) (DELR(J),J=1,NCOL)
c emrl        WRITE(ICNF,*) (DELC(I),I=1,NROW)
c emrl        WRITE(ICNF,*) ((HTOP(J,I),J=1,NCOL),I=1,NROW)
c emrl        WRITE(ICNF,*) (((DZ(J,I,K),J=1,NCOL),I=1,NROW),K=1,NLAY)
c emrl        WRITE(ICNF,*) CINACT,CDRY
c emrl        CLOSE(ICNF)
c emrl      ENDIF
C
C--PROCESS INPUT DATA
C  ==================
C
C--ASSIGN SHARED ICBUND ARRAY TO ALL SPECIES

      DO INDEX=2,NCOMP
        DO K=1,NLAY
          DO I=1,NROW
            DO J=1,NCOL
              ICBUND(J,I,K,INDEX)=ICBUND(J,I,K,1)
            ENDDO
          ENDDO
        ENDDO
      ENDDO
C
C--ASSIGN CINACT TO INACTIVE CELLS AND COPY COLD TO CNEW
      DO INDEX=1,NCOMP
        DO K=1,NLAY
          DO I=1,NROW
            DO J=1,NCOL
              IF(ICBUND(J,I,K,INDEX).EQ.0) COLD(J,I,K,INDEX)=CINACT
              CNEW(J,I,K,INDEX)=COLD(J,I,K,INDEX)
              CADV(J,I,K,INDEX)=COLD(J,I,K,INDEX)
            ENDDO
          ENDDO
        ENDDO
      ENDDO
C
C--SET PRSITY=0 IN CELLS WHERE ICBUND=0
C--AND ENSURE ICBUND=0 IF POROSITY IS ZERO
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(ICBUND(J,I,K,1).EQ.0) THEN
              PRSITY(J,I,K)=0
            ELSEIF(PRSITY(J,I,K).EQ.0) THEN
              DO INDEX=1,NCOMP
                ICBUND(J,I,K,INDEX)=0
              ENDDO
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--CALCULATE COORDINATE ARRAYS XBC, YBC AND ZBC
C--AND MAXIMUN WIDTHS ALONG ROWS, COLUMNS AND LAYERS
      HORIGN=HTOP(1,1)
      XBC(1)=DELR(1)/2.
      YBC(1)=DELC(1)/2.
      DO J=2,NCOL
        XBC(J)=XBC(J-1)+(DELR(J-1)+DELR(J))/2.
      ENDDO
      DO I=2,NROW
        YBC(I)=YBC(I-1)+(DELC(I-1)+DELC(I))/2.
      ENDDO
      XMAX=XBC(NCOL)+DELR(NCOL)/2.
      YMAX=YBC(NROW)+DELC(NROW)/2.
      ZMAX=0
      DO I=1,NROW
        DO J=1,NCOL
          ZBC(J,I,1)=DZ(J,I,1)/2.+(HORIGN-HTOP(J,I))
          ZZ=DZ(J,I,1)
          DO K=2,NLAY
            ZBC(J,I,K)=ZBC(J,I,K-1)+(DZ(J,I,K-1)+DZ(J,I,K))/2.
            ZZ=ZZ+DZ(J,I,K)
          ENDDO
          ZMAX=MAX(ZMAX,ZZ)
        ENDDO
      ENDDO
      WRITE(IOUT,1300) XMAX,YMAX,ZMAX
 1300 FORMAT(/1X,'MAXIMUM LENGTH ALONG THE X (J) AXIS =',G15.7,
     &       /1X,'MAXIMUM LENGTH ALONG THE Y (I) AXIS =',G15.7,
     &       /1X,'MAXIMUM LENGTH ALONG THE Z (K) AXIS =',G15.7)
C
C--INITIALIZE RETARDATION FACTOR ARRAY AND THE MINUMIN
      DO INDEX=1,NCOMP
        DO K=1,NLAY
          DO I=1,NROW
            DO J=1,NCOL
              RETA(J,I,K,INDEX)=1.
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      RFMIN=1.
C
C--RETURN
      RETURN
      END
C
C
      SUBROUTINE BTN4ST(IN,IOUT,NSTP,MXSTP,TSLNGH,DT0,MXSTRN,
     & TTSMULT,TTSMAX,TUNIT)
C *****************************************************************
C THIS SUBROUTINE GETS TIMING INFORMATION FOR EACH STRESS PERIOD.
C *****************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   IN,IOUT,NSTP,MXSTRN,N,MXSTP
      REAL      TSLNGH,DT0,PERLEN,TSMULT,TTSMULT,TTSMAX,EPSILON
      CHARACTER TUNIT*4
      DIMENSION TSLNGH(MXSTP)
      PARAMETER (EPSILON=0.5E-6)
C
C--READ AND PRINT OUT TIMING INFORMATION
      READ(IN,'(F10.0,I10,F10.0)') PERLEN,NSTP,TSMULT
      WRITE(IOUT,22) PERLEN,NSTP,TSMULT
   22 FORMAT(/1X,'LENGTH OF CURRENT STRESS PERIOD =',G15.7,
     & /1X,'NUMBER OF TIME STEPS FOR CURRENT STRESS PERIOD =',I5,
     & /1X,'TIME STEP MULTIPLIER USED IN FLOW SOLUTION =',G15.7)
      IF(NSTP.GT.MXSTP) THEN
        WRITE(*,25)
   25   FORMAT(/1X,'ERROR: MAXIMUM NUMBER OF TIME STEPS EXCEEDED!',
     &   /1X,'INCREASE DIMENSION OF [MXSTP] IN THE MAIN PROGRAM.')
		call stopfile  ! emrl
        STOP
      ENDIF
C
C--IF [TSMULT] IS A NUMBER GREATER THAN ZERO, THE LENGTH OF
C--EACH TIME STEP FOR CURRENT STRESS PERIOD IS CALCULATED BY
C--PROGRAM USING THE GEOMETRIC PROGRESSION.
C--IF TSMULT IS A NUMBER LESS THAN OR EQUAL TO ZERO,
C--READ IN SPECIFIED LENGTH OF EACH TIME STEP.
      IF(TSMULT.LE.0) THEN
        READ(IN,'(8F10.0)') (TSLNGH(N),N=1,NSTP)
        WRITE(IOUT,30) (TSLNGH(N),N=1,NSTP)
   30   FORMAT(1X,'SPECIFIED LENGTH OF EACH TIME STEP:',/1X,8G15.7)
        GOTO 50
      ELSEIF(ABS(TSMULT-1.).LT.ABS(TSMULT+1.)*EPSILON) THEN
        TSLNGH(1)=PERLEN/FLOAT(NSTP)
      ELSE
        TSLNGH(1)=PERLEN*(1.-TSMULT)/(1.-TSMULT**NSTP)
      ENDIF
      DO N=2,NSTP
        TSLNGH(N)=TSLNGH(N-1)*TSMULT
      ENDDO
C
   50 CONTINUE
C
C--READ INITIAL TRANSPORT STEPSIZE AND MAXIMUM NUMBER OF STEPS
      READ(IN,'(F10.0,I10,2F10.0)',ERR=100) DT0,MXSTRN,TTSMULT,TTSMAX
      GOTO 101
  100 BACKSPACE(IN)
      READ(IN,'(F10.0,I10)') DT0,MXSTRN
  101 IF(TTSMULT.LT.1.0) TTSMULT=1.0
      WRITE(IOUT,51) DT0,TUNIT,MXSTRN,TTSMULT,TTSMAX,TUNIT
   51 FORMAT(1X,'USER-SPECIFIED TRANSPORT STEPSIZE =',G15.7,A4
     & /1X,'MAXIMUM NUMBER OF TRANSPORT STEPS ALLOWED ',
     & ' IN ONE FLOW TIME STEP =',I10,
     & /1X,'MULTIPLIER FOR SUCCESSIVE TRANSPORT STEPS ',
     & ' [USED IN IMPLICIT SCHEMES] =',F10.3,
     & /1X,'MAXIMUM TRANSPORT STEP SIZE ',
     & ' [USED IN IMPLICIT SCHEMES] =',G15.7,A4)
C
      IF(DT0.LT.0) WRITE(*,55)
   55 FORMAT(/1X,'NEGATIVE VALUE FOR INPUT VARIABLE [DT0] DETECTED; ',
     & /1X,'MODEL-CALCULATED TRANSPORT STEPSIZE REPLACED WITH [-DT0]',
     & /1X,'REGARDLESS OF STABILITY CONSTRAINTS FOR EXPLICIT SCHEMES',
     & /1X,'OR TRANSPORT STEP MULTIPLIER FOR IMPLICIT SCHEMES.')
C
      RETURN
      END
C
C
      SUBROUTINE BTN4AD(NTRANS,TRNOP,TIME1,TIME2,HT2,DELT,KSTP,NSTP,
     & MXPRS,TIMPRS,DT0,MXSTRN,MIXELM,DTRACK,DTRACK2,
     & PERCEL,DTDISP,DTSSM,DTRCT,RFMIN,NPRS,NPS,DTRANS,PRTOUT,
     & NCOL,NROW,NLAY,NCOMP,ICBUND,CNEW,COLD,CINACT,UPDLHS,IMPSOL,
     & TTSMULT,TTSMAX,KPER,DELR,DELC,DH,PRSITY,SRCONC,RHOB,RETA,
     & PRSITY2,RETA2,ISOTHM,TMASIO,RMASIO,TMASS)
C **********************************************************************
C THIS SUBROUTINE ADVANCES THE TRANSPORT SIMULATION ONE STEP,
C DETERMINING THE STEPSIZE TO BE USED AND WHETHER PRINTOUT IS REQUIRED
C FOR NEXT TRANSPORT STEP. IT ALSO COMPUTES TOTAL MASS IN THE AQUIFER
C AT THE FIRST TRANSPORT STEP OF EACH TRANSPORT LOOP.
C **********************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   NTRANS,KSTP,NSTP,NPS,NPRS,MXPRS,MXSTRN,MIXELM,INDEX,
     &          ICBUND,K,I,J,NCOL,NROW,NLAY,NCOMP,IMPSOL,KPER,ISOTHM
      REAL      TIME1,TIME2,DT0,DTRACK,DTRACK2,PERCEL,DTRANS,DTDISP,
     &          HT2,TIMPRS,DTSSM,DTRCT,DELT,RFMIN,CNEW,COLD,CINACT,
     &          TMASIO,RMASIO,TMASS,DTOLD,TTSMULT,DELR,DELC,DH,
     &          PRSITY,PRSITY2,RETA2,CMML,CMMS,CIML,CIMS,
     &          SRCONC,RHOB,RETA,VOLUME,TTSMAX,EPSILON,TEMP
      DIMENSION ICBUND(NCOL,NROW,NLAY,NCOMP),
     &          CNEW(NCOL,NROW,NLAY,NCOMP),COLD(NCOL,NROW,NLAY,NCOMP),
     &          DELR(NCOL),DELC(NROW),DH(NCOL,NROW,NLAY),
     &          PRSITY(NCOL,NROW,NLAY),RHOB(NCOL,NROW,NLAY),
     &          SRCONC(NCOL,NROW,NLAY,NCOMP),RETA(NCOL,NROW,NLAY,NCOMP),
     &          RETA2(NCOL,NROW,NLAY,NCOMP),PRSITY2(NCOL,NROW,NLAY)
      LOGICAL   TRNOP(10),PRTOUT,UPDLHS
      DIMENSION TIMPRS(MXPRS),TMASIO(122,2,NCOMP),RMASIO(122,2,NCOMP),
     &          TMASS(4,3,NCOMP),TEMP(4)
      PARAMETER (EPSILON=0.5E-6)
C
C--SAVE PREVIOUS TRANSPORT STEPSIZE
      DTOLD=DTRANS
C
C--DETERMINE STEPSIZE FOR NEXT TRANSPORT STEP
      IF(IMPSOL.EQ.0) THEN
        IF(MIXELM.GT.0) THEN
          DTRANS=DELT
          IF(TRNOP(1)) DTRANS=MIN(DTRANS,DTRACK*PERCEL*RFMIN)
          IF(TRNOP(2)) DTRANS=MIN(DTRANS,DTDISP*RFMIN)
          IF(TRNOP(3)) DTRANS=MIN(DTRANS,DTSSM*RFMIN)
          IF(TRNOP(4)) DTRANS=MIN(DTRANS,DTRCT)
        ELSE
          DTRANS=0.
          IF(TRNOP(1)) DTRANS=DTRANS+1./(DTRACK2*PERCEL*RFMIN)
          IF(TRNOP(2)) DTRANS=DTRANS+1./(DTDISP*RFMIN)
          IF(TRNOP(3)) DTRANS=DTRANS+1./(DTSSM*RFMIN)
          IF(TRNOP(4)) DTRANS=DTRANS+1./DTRCT
          DTRANS=1./DTRANS
        ENDIF
        IF(DT0.GT.0.AND.DT0.LT.DTRANS) DTRANS=DT0
      ELSEIF(IMPSOL.EQ.1) THEN
        IF(MIXELM.EQ.0) THEN
          DTRANS=DELT
          IF(TRNOP(1)) DTRANS=MIN(DTRANS,DTRACK*PERCEL*RFMIN)
          IF(DT0.GT.0) DTRANS=DT0
          IF(TTSMULT.GT.1) DTRANS=DTRANS*TTSMULT**(NTRANS-1)
          IF(TTSMAX.GT.0.AND.DTRANS.GT.TTSMAX) DTRANS=TTSMAX
        ELSE
          DTRANS=DELT
          IF(TRNOP(1).AND.MIXELM.GT.0) THEN
            DTRANS=MIN(DTRANS,DTRACK*PERCEL*RFMIN)
          ELSEIF(TRNOP(1).AND.MIXELM.LT.0) THEN
            DTRANS=MIN(DTRANS,DTRACK2*PERCEL*RFMIN)
          ENDIF
          IF(DT0.GT.0.AND.DT0.LT.DTRANS) DTRANS=DT0
          IF(TTSMAX.GT.0.AND.DTRANS.GT.TTSMAX) DTRANS=TTSMAX
        ENDIF
      ENDIF
C
C--IF DT0 NEGATIVE, USE |DT0| AS DEFAULT TRANSPORT STEPSIZE
      IF(DT0.LT.0) DTRANS=ABS(DT0)
C
C--UPDATES TOTAL ELASPED TIME
      TIME1=TIME2
      TIME2=TIME1+DTRANS
C
C--DETERMIN IF PRINTOUT OF SIMULATION RESULTS
C--IS NEEDED FOR NEXT STEP
      PRTOUT=.FALSE.
      IF(NTRANS.EQ.MXSTRN) THEN
        PRTOUT=.TRUE.
      ELSEIF(NPRS.LT.0) THEN
        IF(MOD(NTRANS,-NPRS).EQ.0) PRTOUT=.TRUE.
      ENDIF
C
C--IF TOTAL ELAPSED TIME AT NEXT STEP EXCEEDS TIME AT WHICH
C--PRINTOUT IS REQUESTED, CUT DOWN STEPSIZE
      IF(NPRS.GT.0.AND.NPS.LE.NPRS) THEN
        IF(TIME2.GE.TIMPRS(NPS)) THEN
          IF(TIME2.GT.TIMPRS(NPS)) THEN
            TIME2=TIMPRS(NPS)
            DTRANS=TIME2-TIME1
          ENDIF
          PRTOUT=.TRUE.
          NPS=NPS+1
        ENDIF
      ENDIF
C
C--IF TOTAL ELAPSED TIME AT NEXT TRANSPORT STEP EXCEEDS
C--THE LIMIT OF CURRENT TIME STEP, CUT DOWN STEPSIZE
      IF(TIME2.GE.HT2) THEN
        TIME2=HT2
        DTRANS=TIME2-TIME1
        IF(KSTP.EQ.NSTP) PRTOUT=.TRUE.
      ENDIF
C
      UPDLHS=.TRUE.
      IF(ABS(DTRANS-DTOLD).LT.ABS(DTRANS+DTOLD)*EPSILON
     & .AND. NTRANS.GT.1) UPDLHS=.FALSE.
C
C--PRINT OUT AN IDENTIFYING MSGGAGE
      WRITE(*,70) NTRANS,DTRANS,TIME2
   70 FORMAT(1X,'Transport Step:',I5,3X,'Step Size:',1PG12.4,
     & ' Total Elapsed Time:',G13.5)
C
C--COPY ARRAY [CNEW] TO [COLD]
      DO INDEX=1,NCOMP
        DO K=1,NLAY
          DO I=1,NROW
            DO J=1,NCOL
              IF(ICBUND(J,I,K,INDEX).EQ.0) THEN
                CNEW(J,I,K,INDEX)=CINACT
              ELSE
                COLD(J,I,K,INDEX)=CNEW(J,I,K,INDEX)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDDO
C
C--CLEAR RMASIO ARRAY FOR ACCUMULATING MASS IN/OUT
C--AT NEXT TRANSPORT STEP
      DO INDEX=1,NCOMP
        DO I=1,122
          RMASIO(I,1,INDEX)=0.
          RMASIO(I,2,INDEX)=0.
        ENDDO
      ENDDO
C
C--CALCAULTE TOTAL MASS IN AQUIFER AT THE FIRST TRANSPORT STEP
      IF(NTRANS.GT.1) GOTO 9999
C
C--1: MOBILE-LIQUID   (MML) PHASE
C--2: MOBILE-SORBED   (MMS) PHASE
C--3: IMMOBILE-LIQUID (IML) PHASE
C--4: IMMOBILE-SORBED (IMS) PHASE
C
      DO INDEX=1,NCOMP
        TEMP(1)=0.
        TEMP(2)=0.
        TEMP(3)=0.
        TEMP(4)=0.
        DO K=1,NLAY
          DO I=1,NROW
            DO J=1,NCOL
              IF(ICBUND(J,I,K,INDEX).LE.0) CYCLE
              VOLUME=DELR(J)*DELC(I)*DH(J,I,K)
              CMML=COLD(J,I,K,INDEX)*PRSITY(J,I,K)*VOLUME
              CMMS=0.
              CIML=0.
              CIMS=0.
              IF(ISOTHM.EQ.1) THEN
                CMMS=(RETA(J,I,K,INDEX)-1.)*CMML
              ELSEIF(ISOTHM.GT.1.AND.ISOTHM.LE.4) THEN
                CMMS=SRCONC(J,I,K,INDEX)*RHOB(J,I,K)*VOLUME
              ELSEIF(ISOTHM.GT.4) THEN
                CMMS=(RETA(J,I,K,INDEX)-1.)*CMML
                CIML=PRSITY2(J,I,K)*SRCONC(J,I,K,INDEX)*VOLUME
                CIMS=(RETA2(J,I,K,INDEX)-1.)*CIML
              ENDIF
              TEMP(1)=TEMP(1)+CMML
              TEMP(2)=TEMP(2)+CMMS
              TEMP(3)=TEMP(3)+CIML
              TEMP(4)=TEMP(4)+CIMS
            ENDDO
          ENDDO
        ENDDO
C
C--STORE INITIAL MASS IF THE CURRENT TRANSPORT STEP
C--IS AT THE BEGINNING OF SIMULATION
        IF(KPER*KSTP.EQ.1) THEN
          TMASS(1,1,INDEX)=TEMP(1)
          TMASS(2,1,INDEX)=TEMP(2)
          TMASS(3,1,INDEX)=TEMP(3)
          TMASS(4,1,INDEX)=TEMP(4)
C
C--OTHERWISE DETERMINE AND ACCUMULATE CHANGE IN MASS STORAGE
C--CAUSED BY CHANGE IN SATURATED THICKNESS OF UNCONFINED AQUIFER
        ELSE
          TEMP(1)=TMASS(1,2,INDEX)-TEMP(1)
          TEMP(2)=TMASS(2,2,INDEX)-TEMP(2)
          TEMP(3)=TMASS(3,2,INDEX)-TEMP(3)
          TEMP(4)=TMASS(4,2,INDEX)-TEMP(4)
          TMASS(1,3,INDEX)=TMASS(1,3,INDEX)+TEMP(1)
          TMASS(2,3,INDEX)=TMASS(2,3,INDEX)+TEMP(2)
          TMASS(3,3,INDEX)=TMASS(3,3,INDEX)+TEMP(3)
          TMASS(4,3,INDEX)=TMASS(4,3,INDEX)+TEMP(4)
        ENDIF
      ENDDO
C
 9999 CONTINUE
C
      RETURN
      END
C
C
      SUBROUTINE BTN4SV(NCOL,NROW,NLAY,NCOMP,ICOMP,ICBUND,CNEW,
     & CWGT,CINACT,RMASIO)
C **************************************************************
C THIS SUBROUTINE UPDATES CELL CONCENTRATION AND MASS IN/OUT
C ACCUMULATING ARRAY TO PREPARE FOR SIMULATION AT NEXT STEP.
C **************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   NCOL,NROW,NLAY,NCOMP,ICOMP,J,I,K,ICBUND
      REAL      CNEW,CWGT,CINACT,RMASIO
      DIMENSION ICBUND(NCOL,NROW,NLAY,NCOMP),
     &          CNEW(NCOL,NROW,NLAY,NCOMP),
     &          CWGT(NCOL,NROW,NLAY,NCOMP),RMASIO(122,2,NCOMP)
C
C--COPY CNEW TO CWGT
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(ICBUND(J,I,K,ICOMP).EQ.0) CYCLE
            CWGT(J,I,K,ICOMP)=CNEW(J,I,K,ICOMP)
          ENDDO
        ENDDO
      ENDDO
C
C--CLEAR RMASIO ARRAY FOR ACCUMULATING MASS IN/OUT
C--AT NEXT TRANSPORT STEP
      DO I=1,122
        RMASIO(I,1,ICOMP)=0.
        RMASIO(I,2,ICOMP)=0.
      ENDDO
C
      RETURN
      END
C
C
      SUBROUTINE BTN4BD(KPER,KSTP,NTRANS,NCOL,NROW,NLAY,NCOMP,ICOMP,
     & ISS,ICBUND,DELR,DELC,DH,PRSITY,RETA,CNEW,COLD,RHOB,SRCONC,
     & PRSITY2,RETA2,ISOTHM,DTRANS,TMASIN,TMASOT,ERROR,ERROR2,
     & TMASIO,RMASIO,TMASS)
C **********************************************************************
C THIS SUBROUTINE SUMMARIZES VOLUMETRIC MASS BUDGETS AND CALCULATES
C MASS BALANCE DISCREPANCY SINCE THE BEGINNING OF THE SIMULATION.
C **********************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   KPER,KSTP,NCOL,NROW,NLAY,NCOMP,ICOMP,ICBUND,K,I,J,IQ,
     &          NTRANS,ISOTHM,ISS
      REAL      DELR,DELC,DH,PRSITY,RETA,CNEW,COLD,TMASIN,TMASOT,
     &          TMASIO,RMASIO,DMSTRG,ERROR,ERROR2,TMASS,
     &          SOURCE,SINK,TM1,TM2,DTRANS,RHOB,SRCONC,PRSITY2,RETA2,
     &          CMML,CMMS,CIML,CIMS,VOLUME,STRMAS
      DIMENSION DELR(NCOL),DELC(NROW),DH(NCOL,NROW,NLAY),
     &          RHOB(NCOL,NROW,NLAY),SRCONC(NCOL,NROW,NLAY,NCOMP),
     &          PRSITY(NCOL,NROW,NLAY),RETA(NCOL,NROW,NLAY,NCOMP),
     &          PRSITY2(NCOL,NROW,NLAY),RETA2(NCOL,NROW,NLAY,NCOMP),
     &          CNEW(NCOL,NROW,NLAY,NCOMP),COLD(NCOL,NROW,NLAY,NCOMP),
     &          ICBUND(NCOL,NROW,NLAY,NCOMP),TMASIN(NCOMP),
     &          TMASOT(NCOMP),ERROR(NCOMP),ERROR2(NCOMP),
     &          TMASIO(122,2,NCOMP),RMASIO(122,2,NCOMP),TMASS(4,3,NCOMP)
C
C--CALCULATE SOLUTE AND SORBED MASS STORAGE CHANGES (MOBILE-DOMAIN)
C--FOR THE CURRENT TRANSPORT STEP
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(ICBUND(J,I,K,ICOMP).GT.0.AND.DTRANS.GT.0) THEN
              DMSTRG=(CNEW(J,I,K,ICOMP)-COLD(J,I,K,ICOMP))
     &                   *DELR(J)*DELC(I)*DH(J,I,K)*PRSITY(J,I,K)
              IF(DMSTRG.LT.0) THEN
                RMASIO(119,1,ICOMP)=RMASIO(119,1,ICOMP)-DMSTRG
                RMASIO(120,1,ICOMP)=RMASIO(120,1,ICOMP)
     &           -(RETA(J,I,K,ICOMP)-1.)*DMSTRG
              ELSE
                RMASIO(119,2,ICOMP)=RMASIO(119,2,ICOMP)-DMSTRG
                RMASIO(120,2,ICOMP)=RMASIO(120,2,ICOMP)
     &           -(RETA(J,I,K,ICOMP)-1.)*DMSTRG
              ENDIF
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--ACCUMULATE MASS IN/OUT FOR VARIOUS SINK/SOURCE TERMS AND
C--MASS STOAGE CHANGES SINCE THE BEGINNING OF SIMULATION
      DO IQ=1,122
        TMASIO(IQ,1,ICOMP)=TMASIO(IQ,1,ICOMP)+RMASIO(IQ,1,ICOMP)
        TMASIO(IQ,2,ICOMP)=TMASIO(IQ,2,ICOMP)+RMASIO(IQ,2,ICOMP)
      ENDDO
C
C--DETERMINE TOTAL MASS IN AND OUT
      TMASIN(ICOMP)=0.
      TMASOT(ICOMP)=0.
      DO IQ=1,122
        TMASIN(ICOMP)=TMASIN(ICOMP)+TMASIO(IQ,1,ICOMP)
        TMASOT(ICOMP)=TMASOT(ICOMP)+TMASIO(IQ,2,ICOMP)
      ENDDO
C
C--COMPUTE ACCUMULATIVE DISCREPANCY BETWEEN MASS IN AND OUT
      ERROR(ICOMP)=0.
      IF(ABS(TMASIN(ICOMP))+ABS(TMASOT(ICOMP)).NE.0) THEN
        ERROR(ICOMP)=100.*(TMASIN(ICOMP)+TMASOT(ICOMP))
     &   /(0.5*(ABS(TMASIN(ICOMP))+ABS(TMASOT(ICOMP))))
      ENDIF
C
C--CALCULATE TOTAL MASS IN AQUIFER FOR CURRENT TRANSPORT STEP
      TMASS(1,2,ICOMP)=0.
      TMASS(2,2,ICOMP)=0.
      TMASS(3,2,ICOMP)=0.
      TMASS(4,2,ICOMP)=0.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(ICBUND(J,I,K,ICOMP).LE.0) CYCLE
            VOLUME=DELR(J)*DELC(I)*DH(J,I,K)
            CMML=CNEW(J,I,K,ICOMP)*PRSITY(J,I,K)*VOLUME
            CMMS=0.
            CIML=0.
            CIMS=0.
            IF(ISOTHM.EQ.1) THEN
              CMMS=(RETA(J,I,K,ICOMP)-1.)*CMML
            ELSEIF(ISOTHM.GT.1.AND.ISOTHM.LE.4) THEN
              CMMS=SRCONC(J,I,K,ICOMP)*RHOB(J,I,K)*VOLUME
            ELSEIF(ISOTHM.GT.4) THEN
              CMMS=(RETA(J,I,K,ICOMP)-1.)*CMML
              CIML=PRSITY2(J,I,K)*SRCONC(J,I,K,ICOMP)*VOLUME
              CIMS=(RETA2(J,I,K,ICOMP)-1.)*CIML
            ENDIF
            TMASS(1,2,ICOMP)=TMASS(1,2,ICOMP)+CMML
            TMASS(2,2,ICOMP)=TMASS(2,2,ICOMP)+CMMS
            TMASS(3,2,ICOMP)=TMASS(3,2,ICOMP)+CIML
            TMASS(4,2,ICOMP)=TMASS(4,2,ICOMP)+CIMS
          ENDDO
        ENDDO
      ENDDO
C
C--COMPUTE TOTAL SOURCE AND SINK EXCLUDING MASS
C--FROM OR INTO FLUID-STORAGE IN TRANSIENT FLOW FIELD
      SOURCE=0.
      SINK=0.
      DO IQ=1,117
        SOURCE=SOURCE+TMASIO(IQ,1,ICOMP)
        SINK=SINK+TMASIO(IQ,2,ICOMP)
      ENDDO
C
C--GET SUM OF TOTAL SOURCE AND INITIAL MASS
      TM1=ABS(SOURCE)+TMASS(1,1,ICOMP)+TMASS(2,1,ICOMP)
     &   +TMASS(3,1,ICOMP)+TMASS(4,1,ICOMP)
C
C--GET SUM OF TOTAL SINK AND CURRENT MASS
      TM2=ABS(SINK)+TMASS(1,2,ICOMP)+TMASS(2,2,ICOMP)
     &   +TMASS(3,2,ICOMP)+TMASS(4,2,ICOMP)
C
C--CORRECT FOR NET MASS FROM/INTO FLUID-STORAGE
      STRMAS=(TMASIO(118,1,ICOMP)+TMASIO(118,2,ICOMP))
     &     -(TMASS(1,3,ICOMP)+TMASS(2,3,ICOMP)
     &      +TMASS(3,3,ICOMP)+TMASS(4,3,ICOMP))
      TM1=TM1+STRMAS
C
C--COMPUTE ALTERNATVE MEASURE OF MASS DISCREPANCY
      ERROR2(ICOMP)=0.
      IF( TM1+TM2 .NE.0. ) THEN
        ERROR2(ICOMP)=100.*(TM1-TM2)/(0.5*(TM1+TM2))
      ENDIF
C
C--RETURN
      RETURN
      END
C
C
C emrl - added ispc to subroutine below
      SUBROUTINE BTN4OT(NCOL,NROW,NLAY,KPER,KSTP,NTRANS,NCOMP,ICOMP,
     & IOUT,IOBS,IUCN,ispc,IMAS,ICBM,MXOBS,NOBS,NPROBS,LOCOBS,ICBUND,
     & TIME2,CNEW,MIXELM,NCOUNT,NPINS,NRC,NPCHEK,ISOTHM,RETA,TMASIN,
     & TMASOT,ERROR,ERROR2,TRNOP,TUNIT,MUNIT,PRTOUT,TMASIO,RMASIO,TMASS)
C **********************************************************************
C THIS SUBROUTINE SAVES SIMULATION RESULTS IN THE STANDARD OUTPUT FILE
C AND VARIOUS OPTIONAL OUTPUT FILES, ACCORDING TO THE OUTPUT CONTROL
C OPTIONS SPECIFIED IN THE BASIC TRANSPORT INPUT FILE.
C **********************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   NCOL,NROW,NLAY,KPER,KSTP,NTRANS,IOUT,IOBS,IUCN,IMAS,
     &          ICBM,MIXELM,NCOUNT,NPINS,NRC,ISOTHM,NPCHEK,MXOBS,NOBS,
     &          NPROBS,LOCOBS,IFMTCN,IFMTNP,IFMTRF,IFMTDP,K,I,J,N,
     &          NPRMAS,ICBUND,NCOMP,ICOMP
      REAL      TIME2,CNEW,RETA,TMASIN,TMASOT,ERROR,TMASIO,RMASIO,
     &          TMASS,ERROR2,SOURCE,SINK,STRMAS,TOTMAS
      LOGICAL   TRNOP(10),FWEL,FDRN,FRCH,FEVT,FRIV,FGHB,
     &          FSTR,FRES,FFHB,FIBS,FTLK,FLAK,FMAW,FDRT,FETS,FUSR(3),
     &          PRTOUT,CHKMAS,SAVUCN,SAVCBM
      CHARACTER TEXT*16,TUNIT*4,MUNIT*4
      DIMENSION CNEW(NCOL,NROW,NLAY,NCOMP),NPCHEK(NCOL,NROW,NLAY,NCOMP),
     &          RETA(NCOL,NROW,NLAY,NCOMP),
     &          ICBUND(NCOL,NROW,NLAY,NCOMP),LOCOBS(3,MXOBS),
     &          TMASIO(122,2,NCOMP),RMASIO(122,2,NCOMP),
     &          TMASS(4,3,NCOMP),NCOUNT(NCOMP),NPINS(NCOMP),NRC(NCOMP),
     &          TMASIN(NCOMP),TMASOT(NCOMP),ERROR(NCOMP),ERROR2(NCOMP)
      COMMON /FC/FWEL,FDRN,FRCH,FEVT,FRIV,FGHB,
     &          FSTR,FRES,FFHB,FIBS,FTLK,FLAK,FMAW,FDRT,FETS,FUSR
      COMMON /OC/IFMTCN,IFMTNP,IFMTRF,IFMTDP,SAVUCN,SAVCBM,CHKMAS,NPRMAS
C-------emrl start
      integer ITS,ISTATE,ISPC
C-------emrl end
C
C--PRINT OUT CONCENTRATIONS AT SPECIFIED OBSERVATION POINTS
C--TO FILE [MT3Dnnn.OBS] IF REQUESTED.
      IF(NOBS.GT.0.AND.(MOD(NTRANS-1,NPROBS).EQ.0.OR.PRTOUT)) THEN
        IF(NOBS.LE.16) THEN
          WRITE(IOBS+ICOMP,1000) NTRANS,TIME2,(CNEW(LOCOBS(3,N),
     &     LOCOBS(2,N),LOCOBS(1,N),ICOMP),N=1,NOBS)
        ELSE
          WRITE(IOBS+ICOMP,1010) NTRANS,TIME2,(CNEW(LOCOBS(3,N),
     &     LOCOBS(2,N),LOCOBS(1,N),ICOMP),N=1,NOBS)
        ENDIF
      ENDIF
 1000 FORMAT(1X,I5,1X,1PG13.5,1X,16(G13.5,1X))
 1010 FORMAT(1X,I5,1X,1PG13.5,1X,16(G13.5,1X)/(1X,20X,16(G13.5,1X)))
C
C--WRITE A ONE-LINE SUMMARY OF MASS BALANCE
C--TO FILE [MT3Dnnn.MAS] IF REQUESTED
      IF(CHKMAS.AND.(MOD(NTRANS-1,NPRMAS).EQ.0.OR.PRTOUT)) THEN
        SOURCE=0.
        SINK=0.
        DO N=1,117
          SOURCE=SOURCE+TMASIO(N,1,ICOMP)
          SINK=SINK+TMASIO(N,2,ICOMP)
        ENDDO
        STRMAS=(TMASIO(118,1,ICOMP)+TMASIO(118,2,ICOMP))
     &        -(TMASS(1,3,ICOMP)+TMASS(2,3,ICOMP)
     &        +TMASS(3,3,ICOMP)+TMASS(4,3,ICOMP))
        TOTMAS=TMASS(1,2,ICOMP)+TMASS(2,2,ICOMP)
     &        +TMASS(3,2,ICOMP)+TMASS(4,2,ICOMP)
        WRITE(IMAS+ICOMP,1012) TIME2,TMASIN(ICOMP),TMASOT(ICOMP),
     &   SOURCE,SINK,STRMAS,TOTMAS,
     &   ERROR(ICOMP),ERROR2(ICOMP)
      ENDIF
 1012 FORMAT(1X,1P,7(G13.5,1X),4X,G10.3,5X,G10.3)
C
C--SAVE CELL CONCENTRATIONS TO UNFORMATTED FILE [MT3Dnnn.UCN]
C--IF REQUESTED
c emrl      IF(SAVUCN .AND. PRTOUT) THEN
c emrl        TEXT='CONCENTRATION'
c emrl        DO K=1,NLAY
c emrl          WRITE(IUCN+ICOMP) NTRANS,KSTP,KPER,TIME2,TEXT,NCOL,NROW,K
c emrl          WRITE(IUCN+ICOMP) ((CNEW(J,I,K,ICOMP),J=1,NCOL),I=1,NROW)
c emrl        ENDDO
c emrl      ENDIF
C---- emrl start
      IF(SAVUCN .AND. PRTOUT) THEN
        ITS = 200
        ISTATE = 1
        WRITE(ISPC+ICOMP) ITS,ISTATE,TIME2
        WRITE(ISPC+ICOMP) ((((ABS(ICBUND(J,I,K,ICOMP))/
     &    MAX(ABS(ICBUND(J,I,K,ICOMP)),1)),J=1,NCOL),I=1,NROW),K=1,NLAY)
        WRITE(ISPC+ICOMP) (((CNEW(J,I,K,ICOMP),J=1,NCOL),
     &    I=1,NROW),K=1,NLAY)
      ENDIF
C---- emrl end
C
C--WRITE SIMULATION RESULTS AND MASS BUDGET TERMS
C--TO THE STANDARD OUTFILE IF NEEDED
      IF(.NOT.PRTOUT) GOTO 9999
C
C--PRINT A HEADER
      WRITE(IOUT,1019) ICOMP
 1019 FORMAT(///1X,55('>'),'FOR COMPONENT NO.',I3.2,55('<'))
C
      WRITE(IOUT,1020) NTRANS,TIME2,TUNIT
 1020 FORMAT(//44X,43('-')/54X,'TRANSPORT STEP NO.',I5/44X,43('-')
     & //1X,'TOTAL ELAPSED TIME SINCE BEGINNING OF SIMULATION =',
     & G15.7,A4/1X,69('.'))
C
C--PRINT CELL CONCENTRATIONS IF NEEDED
      IF(IFMTCN.EQ.0) GOTO 40
C
      TEXT='CONCENTRATIONS'
      DO K=1,NLAY
        CALL RPRINT(CNEW(1,1,K,ICOMP),TEXT,
     &   NTRANS,KSTP,KPER,NCOL,NROW,K,IFMTCN,IOUT)
      ENDDO
C
C--PRINT NONLINEAR RETARDATION FACTOR IF NEEDED
   40 IF(.NOT.TRNOP(4)) GOTO 50
      IF(ISOTHM.NE.2.AND.ISOTHM.NE.3) GOTO 50
      IF(IFMTRF.EQ.0) GOTO 50
C
      TEXT='RETARD. FACTOR'
      DO K=1,NLAY
        CALL RPRINT(RETA(1,1,K,ICOMP),TEXT,
     &   NTRANS,KSTP,KPER,NCOL,NROW,K,IFMTRF,IOUT)
      ENDDO
C
C--PRINT PARTICLE USAGE INFORMATION IF TRNOP(1)=T AND MIXELM>0
   50 IF(.NOT.TRNOP(1)) GOTO 70
      IF(MIXELM.LE.0) GOTO 70
C
      WRITE(IOUT,1030) NCOUNT(ICOMP),NPINS(ICOMP),NRC(ICOMP)
 1030 FORMAT(/1X,'TOTAL PARTICLES USED IN THE CURRENT STEP =',I10
     &       /1X,'PARTICLES ADDED AT BEGINNING OF THE STEP =',I10
     &       /1X,'PARTICLES REMOVED AT END OF LAST STEP    =',I10)
C
C--PRINT PARTICLE NUMBER PER CELL IF NEEDED
      IF(IFMTNP.EQ.0) GOTO 70
C
      TEXT='PARTICLE NUMBER '
      DO K=1,NLAY
        CALL IPRINT(NPCHEK(1,1,K,ICOMP),TEXT,
     &   NTRANS,KSTP,KPER,NCOL,NROW,K,IFMTNP,IOUT)
      ENDDO
C
C--PRINT OUT ACCUMULATIVE MASS BALANCE INFORMATION
   70 WRITE(IOUT,1110) NTRANS,KSTP,KPER
      WRITE(IOUT,1114)
      WRITE(IOUT,1122) TMASIO(6,1,ICOMP),TMASIO(6,2,ICOMP)
      WRITE(IOUT,1120) TMASIO(1,1,ICOMP),TMASIO(1,2,ICOMP)
      IF(FWEL) WRITE(IOUT,1130) TMASIO(2,1,ICOMP),TMASIO(2,2,ICOMP)
      IF(FDRN) WRITE(IOUT,1140) TMASIO(3,1,ICOMP),TMASIO(3,2,ICOMP)
      IF(FRIV) WRITE(IOUT,1150) TMASIO(4,1,ICOMP),TMASIO(4,2,ICOMP)
      IF(FGHB) WRITE(IOUT,1160) TMASIO(5,1,ICOMP),TMASIO(5,2,ICOMP)
      IF(FRCH) WRITE(IOUT,1162) TMASIO(7,1,ICOMP),TMASIO(7,2,ICOMP)
      IF(FEVT) WRITE(IOUT,1164) TMASIO(8,1,ICOMP),TMASIO(8,2,ICOMP)
      IF(TMASIO(15,1,ICOMP)-TMASIO(15,2,ICOMP).NE.0)
     &   WRITE(IOUT,1165) TMASIO(15,1,ICOMP),TMASIO(15,2,ICOMP)
C
      IF(FSTR) WRITE(IOUT,2100) TMASIO(21,1,ICOMP),TMASIO(21,2,ICOMP)
      IF(FRES) WRITE(IOUT,2102) TMASIO(22,1,ICOMP),TMASIO(22,2,ICOMP)
      IF(FFHB) WRITE(IOUT,2104) TMASIO(23,1,ICOMP),TMASIO(23,2,ICOMP)
      IF(FIBS) WRITE(IOUT,2106) TMASIO(24,1,ICOMP),TMASIO(24,2,ICOMP)
      IF(FTLK) WRITE(IOUT,2108) TMASIO(25,1,ICOMP),TMASIO(25,2,ICOMP)
      IF(FLAK) WRITE(IOUT,2110) TMASIO(26,1,ICOMP),TMASIO(26,2,ICOMP)
      IF(FMAW) WRITE(IOUT,2112) TMASIO(27,1,ICOMP),TMASIO(27,2,ICOMP)
      IF(FDRT) WRITE(IOUT,2114) TMASIO(28,1,ICOMP),TMASIO(28,2,ICOMP)
      IF(FETS) WRITE(IOUT,2116) TMASIO(29,1,ICOMP),TMASIO(29,2,ICOMP)
      IF(FUSR(1)) WRITE(IOUT,2200) TMASIO(51,1,ICOMP),TMASIO(51,2,ICOMP)
      IF(FUSR(2)) WRITE(IOUT,2202) TMASIO(52,1,ICOMP),TMASIO(52,2,ICOMP)
      IF(FUSR(3)) WRITE(IOUT,2204) TMASIO(53,1,ICOMP),TMASIO(53,2,ICOMP)
C
      IF(TRNOP(4)) WRITE(IOUT,1166) TMASIO(9,1,ICOMP),TMASIO(9,2,ICOMP)
      WRITE(IOUT,1170) TMASIO(119,1,ICOMP),TMASIO(119,2,ICOMP)
      IF(TRNOP(4).AND.ISOTHM.GT.0)
     &   WRITE(IOUT,1172) TMASIO(120,1,ICOMP),TMASIO(120,2,ICOMP)
      IF(TRNOP(4).AND.ISOTHM.GT.4) THEN
        WRITE(IOUT,1173)
        WRITE(IOUT,1174) TMASIO(10,1,ICOMP),TMASIO(10,2,ICOMP)
        WRITE(IOUT,1175) TMASIO(121,1,ICOMP),TMASIO(121,2,ICOMP)
        WRITE(IOUT,1176) TMASIO(122,1,ICOMP),TMASIO(122,2,ICOMP)
      ENDIF
      WRITE(IOUT,1180) TMASIN(ICOMP),MUNIT,TMASOT(ICOMP),MUNIT,
     &   TMASIN(ICOMP)+TMASOT(ICOMP),ERROR(ICOMP)
C
 1110 FORMAT(/21X,'CUMMULATIVE MASS BUDGETS AT END OF TRANSPORT STEP',
     & I5,', TIME STEP',I5,', STRESS PERIOD',I5/21X,90('-'))
 1114 FORMAT(/30X,24X,7X,'IN',8X,13X,6X,'OUT',
     &      /30X,24X,16('-'),13X,16('-'))
 1122 FORMAT(30X,' CONSTANT CONCENTRATION: ',G15.7,13X,G15.7)
 1120 FORMAT(30X,'          CONSTANT HEAD: ',G15.7,13X,G15.7)
 1130 FORMAT(30X,'                  WELLS: ',G15.7,13X,G15.7)
 1140 FORMAT(30X,'                 DRAINS: ',G15.7,13X,G15.7)
 1150 FORMAT(30X,'                 RIVERS: ',G15.7,13X,G15.7)
 1160 FORMAT(30X,'HEAD-DEPENDENT BOUNDARY: ',G15.7,13X,G15.7)
 1162 FORMAT(30X,'               RECHARGE: ',G15.7,13X,G15.7)
 1164 FORMAT(30X,'     EVAPOTRANSPIRATION: ',G15.7,13X,G15.7)
 1165 FORMAT(30X,'           MASS LOADING: ',G15.7,13X,G15.7)
C
 2100 FORMAT(30X,'     STREAMFLOW ROUTING: ',G15.7,13X,G15.7)
 2102 FORMAT(30X,'              RESERVOIR: ',G15.7,13X,G15.7)
 2104 FORMAT(30X,' FLOW AND HEAD BOUNDARY: ',G15.7,13X,G15.7)
 2106 FORMAT(30X,'  INTERBED STORAGE FLOW: ',G15.7,13X,G15.7)
 2108 FORMAT(30X,'     TRANSIENT LEACKAGE: ',G15.7,13X,G15.7)
 2110 FORMAT(30X,'                   LAKE: ',G15.7,13X,G15.7)
 2112 FORMAT(30X,'     MULTI-AQUIFER WELL: ',G15.7,13X,G15.7)
 2114 FORMAT(30X,' DRAIN WITH RETURN FLOW: ',G15.7,13X,G15.7)
 2116 FORMAT(30X,'           SEGMENTED ET: ',G15.7,13X,G15.7)
 2200 FORMAT(30X,'      USER-DEFINED NO.1: ',G15.7,13X,G15.7)
 2202 FORMAT(30X,'      USER-DEFINED NO.2: ',G15.7,13X,G15.7)
 2204 FORMAT(30X,'      USER-DEFINED NO.3: ',G15.7,13X,G15.7)
C
 1166 FORMAT(30X,'DECAY OR BIODEGRADATION: ',G15.7,13X,G15.7)
 1170 FORMAT(30X,'  MASS STORAGE (SOLUTE): ',G15.7,13X,G15.7)
 1172 FORMAT(30X,'  MASS STORAGE (SORBED): ',G15.7,13X,G15.7)
 1173 FORMAT(30X,'....immobile domain....')
 1174 FORMAT(30X,'DECAY OR BIODEGRADATION: ',G15.7,13X,G15.7)
 1175 FORMAT(30X,'  MASS STORAGE (SOLUTE): ',G15.7,13X,G15.7)
 1176 FORMAT(30X,'  MASS STORAGE (SORBED): ',G15.7,13X,G15.7)
 1180 FORMAT(28X,75('-'),
     &      /30X,'                [TOTAL]: ',G15.7,1X,A4,8X,G15.7,
     & 1X,A4//40X,'                  NET (IN - OUT): ',G15.7,
     &       /40X,'           DISCREPANCY (PERCENT): ',G15.7)
C
C--RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE BTN4FM(NCOL,NROW,NLAY,NCOMP,ICOMP,ICBUND,CADV,COLD,
     & RETA,PRSITY,DELR,DELC,DZ,DTRANS,A,RHS,NODES,UPDLHS,NCRS,MIXELM)
C *********************************************************************
C THIS SUBROUTINE INITIALIZES ALL MATRICES FOR THE IMPLICIT SCHEME.
C *********************************************************************
C last modified: 08-12-2001
C
      IMPLICIT  NONE
      INTEGER   NCOL,NROW,NLAY,NCOMP,ICOMP,J,I,K,ICBUND,NODES,N,NRC,
     &          NCRS,NSIZE,L,MIXELM
      REAL      CADV,COLD,A,RHS,PRSITY,DTRANS,DZ,RETA,DELR,DELC,TEMP
      LOGICAL   UPDLHS
      DIMENSION ICBUND(NODES,NCOMP),CADV(NODES,NCOMP),COLD(NODES,NCOMP),
     &          RETA(NODES,NCOMP),PRSITY(NODES),DZ(NODES),DELR(NCOL),
     &          DELC(NROW),A(19*NODES),RHS(NODES)
      COMMON   /GCGIDX/L(19)
C
C--GET RIGHT-HAND-SIDE ARRAY [RHS]
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            N=(K-1)*NCOL*NROW + (I-1)*NCOL + J
            IF(MIXELM.EQ.0) THEN
              TEMP=COLD(N,ICOMP)
            ELSE
              TEMP=CADV(N,ICOMP)
            ENDIF
            IF(ICBUND(N,ICOMP).LE.0) THEN
              RHS(N)=-TEMP
            ELSE
              RHS(N)=-TEMP*RETA(N,ICOMP)/DTRANS*PRSITY(N)
     &               *DELR(J)*DELC(I)*DZ(N)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RETURN IF COEFF MATRIX [A] NOT TO BE UPDATED
      IF(.NOT.UPDLHS) GOTO 999
C
C--RESET COEFF MATRIX [A]
C--IF ALL CROSS TERMS ARE INVOLVED, ARRAY [A] HAS
C--LENGTH OF 19 * NODES
      IF(NCRS.EQ.1) THEN
         NSIZE=19*NODES
C
C--OTHERWISE IT HAS LENGTH OF 7 * NODES
      ELSE
         NSIZE=7*NODES
      ENDIF
C
C--CLEAR THE ARRAY
      DO I=1,NSIZE
        A(I)=0.
      ENDDO
C
C--LOOP THROUGH ALL CELLS AND RESET A
      N=0
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            N=N+1
C
C--IF INACTIVE OR CONSTANT CELL
            IF(ICBUND(N,ICOMP).LE.0) THEN
               A(N)=-1.
            ELSE
               A(N)=-RETA(N,ICOMP)/DTRANS*PRSITY(N)
     &              *DELR(J)*DELC(I)*DZ(N)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--CALCULATE MATRIX INDICES FOR THE GCG SOLVER
      NRC = NROW*NCOL
      L(1) =  0
      L(2) = -NRC
      L(3) =  NRC
      L(4) = -NCOL
      L(5) =  NCOL
      L(6) = -1
      L(7) =  1
      L(8) = -NCOL-NRC
      L(9) = -1-NRC
      L(10)= 1-NRC
      L(11)= NCOL-NRC
      L(12)=-NCOL+NRC
      L(13)=-1+NRC
      L(14)= 1+NRC
      L(15)= NCOL+NRC
      L(16)=-1-NCOL
      L(17)= 1-NCOL
      L(18)=-1+NCOL
      L(19)= 1+NCOL
C
C--NORMAL RETURN
  999 RETURN
      END
C
C
      LOGICAL FUNCTION UNIFOR(A,NC,NR,NL)
C ***************************************************
C THIS FUNCTION CHECKS WHETHER ELEMENTS IN AN ARRAY
C ARE UNIFORM.
C ***************************************************
C last modified: 06-23-1998
C
      IMPLICIT  NONE
      INTEGER   NC,NR,NL,J,I,K
      REAL      A,AI,EPSILON
      PARAMETER (EPSILON=0.5E-6)
      DIMENSION A(NC,NR,NL)
C
C--GET THE 1ST ELEMENT
      AI=A(1,1,1)
C
C--COMPARE REST OF ELEMENTS WITH THE 1ST ELEMENT
      DO K=1,NL
        DO I=1,NR
          DO J=1,NC
            IF(ABS(A(J,I,K)-AI).GT.ABS(A(J,I,K)+AI)*EPSILON) THEN
              UNIFOR=.FALSE.
              RETURN
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--IF ALL ELEMENTS ARE EQUAL, SET [UNIFOR] TO T
      UNIFOR=.TRUE.
      RETURN
      END
